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Abstract: The hard-thermal- loop perturbation theory (HTLpt) framework is used to 
calculate the thermodynamic functions of a quark- gluon plasma to three-loop order. This 
is the highest order accessible by finite temperature perturbation theory applied to a non- 
Abelian gauge theory before the high-temperature infrared catastrophe. All ultraviolet 
divergences are eliminated by renormalization of the vacuum, the HTL mass parameters, 
and the strong coupling constant. After choosing a prescription for the mass parameters, 
the three-loop results for the pressure and trace anomaly are found to be in very good 
agreement with recent lattice data down to T ~ 2 — 3 Tc, which are temperatures accessible 
by current and forthcoming heavy- ion collision experiments. 
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1 Introduction 

The last generation of relativistic heavy-ion collision experiments exceeded the energy den- 
sity necessary for the formation of a quark-gluon plasma, motivating the development of a 
quantitative framework with which to calculate the properties of this new state of matter. 
At Brookhaven National Laboratory (RHIC), New York, USA, the initial temperatures 
were up to twice the critical temperature for deconfinement,^, Tc ~ 170 MeV. This trans- 
lates to a strong coupling constant of as(/i = 27r x 170 MeV) = 5'^/(47r) ~ 0.43 where ^ 
is the renormalization scale and it relates to the temperature by = 27rT. The upcoming 
experiments at CERN (LHC), Geneva, Switzerland, are expected to yield initial tempera- 
tures of 4 — 6 Tc, driving the running coupling further down. Initially, the hope was that the 
asymptotic freedom of QCD would allow calculations using a perturbative expansion in the 
coupling. Utilizing a weak-coupling expansion in the coupling constant g to calculate the 
thermodynamic functions of QCD has a long history [2-9], and the pressure is now known 
through order g^logg for non-Abelian gauge theories [10]. Unfortunately, for all but tiny 
coupling constants, and thus astronomically high temperatures, these expansions converge 

^Note that the deconfinement transition is actually an analytic crossover [1] and Tc represents a tem- 
perature around which the thermodynamic quantities change quickly. 
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Figure 1. Weak-coupling expansion for the scaled QCD pressure with Nf — 3. Shaded bands show 
the result of varying the renormalization scale /i by a factor of 2 around /i = 2ttT. 

very poorly, and show large dependence on the renormalization scale. In figure 1, we show 
the weak-coupling expansion for the QCD pressure with Nj = 3 normalized to that of an 
ideal gas through order g^. The various approximations oscillate wildly and show no signs 
of convergence in the temperature range shown. The bands are obtained by varying the 
renormalization scale by a factor of 2 around the value = 2itT and we use three-loop 
running for Og [11] with Ajj^{Nf = 3) = 344 MeV [12]. This oscillating behavior is generic 
for hot field theories, and not specific to QCD. The instability is thought to be caused 
by plasma effects such as screening of electric fields and Landau damping. This calls for 
a nonperturbative approach, or a reorganization of the perturbative expansion that takes 
such effects into account. Furthermore, data from RHIC suggest that the matter created 
behaves more like a strongly coupled fluid with small viscosity [13-17], inspiring the de- 
velopment of strongly coupled formalisms, perhaps the most successful being those based 
on the Anti-de-Sitter space/conformal field theory (AdS/CFT) correspondence proposed 
by Maldacena [18]. However, work on the perturbative side has shown that observables 
like jet quenching [19, 20] and elliptic flow [21] can also be described by a perturbative 
setup. Looking forward to the upcoming heavy-ion experiments scheduled to take place 
at the LHC, it is therefore important to know if, at the higher temperatures generated, 
one should expect a strongly-coupled (liquid) or weakly-coupled (plasma) description to 
be more appropriate. The key question is, will the generated matter behave more like a 
plasma of quasiparticles at these higher temperatures. 

Of course, another approach is lattice gauge theory, which is a nonperturbative dis- 
crete space-time framework that is the closest one can currently get to a first-principles 
calculation with realistic parameters. However, the Monte Carlo methods used currently 
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restrict lattice gauge theory to the study of systems with low chemical potential, and re- 
liable calculations of dynamical quantities have proved difficult. A quark-gluon plasma 
created in a heavy-ion collision will be out of thermal equilibrium and will have nonzero 
baryon density, making these restrictions quite severe. Still, the part of the phase diagram 
that is accessible by lattice gauge theory becomes a clean testing ground for the various 
approaches to analytical QCD calculations. In this paper we will be comparing our new 
"reorganized" perturbative results for the thermodynamic functions to recent lattice data 
from refs. [22] and [23]. 

There are several ways of systematically reorganizing the perturbative expansion [24- 
26]. In screened perturbation theory (SPT) [27-31] which was inspired by variational 
perturbation theory (VPT) [32-37], a mass term is added to the free Lagrangian, and 
then subtracted as an interaction term. Unfortunately, this reorganization cannot be di- 
rectly applied to gauge theories, since the added local mass term would violate gauge 
invariance [38-40]. Instead, one can add and subtract a hard-thermal-loop (HTL) effective 
action [41] which modifies the propagators and vertices systematically and self-consistently, 
so that the reorganization is manifestly gauge invariant [42-44] . The mass parameters m^i 
and ruq are introduced, and identified with the Debye screening mass and the thermal 
quark mass at leading order, respectively, in order to reproduce the high temperature limit 
of thermal QCD; however, at higher orders it is necessary to introduce a prescription for 
the subleading contributions to the mass parameters. The resulting HTL perturbation 
theory (HTLpt) framework can be applied to static as well as dynamic quantities. 

The HTLpt framework has recently been applied to QED [45] , where the calculation of 
the thermodynamic functions is carried out through next-to-next-to- leading-order (NNLO). 
The thermodynamic functions of QCD have been calculated to next-to-leading-order (NLO) 
[46, 47]. For pure-glue QCD the thermodynamic functions were recently calculated to 
NNLO in [48, 49]. Our paper builds upon the results of NNLO QED and pure-glue QCD, 
and we now include NNLO contributions from quark and quark-gluon interaction diagrams 
to express the thermodynamic functions of full QCD to NNLO. The calculation is presented 
in some detail. For a short letter with just the final results, we refer to ref. [50]. Our results 
indicate that the lattice data at temperatures T ~ 2 Tc are consistent with the quasiparticle 
picture. This is a nontrivial result since, in this temperature regime, the QCD coupling 
constant is neither infinitesimally weak nor infinitely strong with (7 ~ 2, or equivalently 
Us ~ 0.3. Therefore, we have a crucial test of the quasiparticle picture in the intermediate 
coupling regime. 

After the introduction, we begin in section 2 with a summary of HTLpt applied to 
QCD. Section 3 contains expressions for the diagrams contributing to the thermodynamic 
potential up to NNLO. In section 4 these diagrams are reduced to scalar sum- integrals and 
expanded in powers of uid/T and niq/T, keeping all terms that contribute through order 
if the mass parameters are taken to be of order g at leading order (LO). In section 5 
the calculated diagrams are gathered to obtain the renormalized NNLO thermodynamic 
potential. The results are then presented in section 6, and compared to recent lattice data. 
Due to difficulties with the normal variational approach, there is special emphasis placed 
on the determination of the mass parameters m^) and niq. In section 7 we summarize. 
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2 HTL perturbation theory 



The Lagrangian density for QCD in Minkowski space is 

Cqcd = -^Tr + ii^l^D^i; + £gf + £gh + A£qcd , (2.1) 

where the gluon field strength is G'^'^ = d'^A'^ — d'^A'^ — iglA^^, A'^], the term with the quark 
fields ip contains an implicit sum over the Nf quark flavors, and the covariant derivative is 
jjfi — Qfj- _ igAfJ-, The ghost term £gh depends on the gauge- fixing term £gf . In this paper 
we choose the class of covariant gauges where the gauge-fixing term is 



^gf = -^Tr 



{df^An^ . (2.2) 



The perturbative expansion in powers of g generates ultraviolet divergences. The 
renormalizability of perturbative QCD guarantees that all divergences in physical quan- 
tities can be removed by renormalization of the coupling constant Ug = (7^/(47r) and the 
necessary counterterms are represented by A£qcd in the Lagrangian (2.1). There is no 
need for wavefunction renormalization, because physical quantities are independent of the 
normalization of the field. There is also no need for renormalization of the gauge parameter, 
because physical quantities are independent of the gauge parameter. 

Hard-thermal-loop perturbation theory (HTLpt) is a reorganization of the perturbation 
series for thermal QCD. The Lagrangian density is written as 

C = (£qcd + /:htl) ^ + A£htl. (2.3) 

The HTL improvement term is 

W = - ^(1 - ^)"^DTr {g,^ (jP^) + (1 - (^) ^' ^^-^^ 

where = (1, y) is a light-like four-vector, and (. . .)y represents the average over the direc- 
tions of y. The term (2.4) has the form of the effective Lagrangian that would be induced 
by a rotationally invariant ensemble of charged sources with infinitely high momentum and 
modifies the propagators and vertices self-consistently so that the reorganization is man- 
ifestly gauge invariant [38, 51-55]. The parameter mo can be identified with the Debye 
screening mass, and rriq with the thermal quark mass to account for the screening effects. 
HTLpt is defined by treating 5 as a formal expansion parameter. By coupling the HTL 
improvement term (2.4) to the QCD Lagrangian (2.1), HTLpt systematically shifts the 
perturbative expansion from being around an ideal gas of massless particles which is the 
physical picture of the weak-coupling expansion, to being around a gas of massive quasi- 
particles which are the more appropriate physical degrees of freedom at high temperature. 

Physical observables are calculated in HTLpt by expanding them in powers of 5, trun- 
cating at some specified order, and then setting 5 = 1. This defines a reorganization of the 
perturbation series in which the effects of and rrig terms in (2.4) are included to all 
orders but then systematically subtracted out at higher orders in perturbation theory by 
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the 6m'j-) and 6mg terms in (2.4). If we set (5 = 1, the HTLpt Lagrangian (2.3) reduces to 
the QCD Lagrangian (2.1). 

If the expansion in 6 could be calculated to all orders the final result would not depend 
on niD and niq when we set 5=1. However, any truncation of the expansion in 5 produces 
results that depend on m/j and mg. Some prescription is required to determine mo and 
ruq as a function of T and as- We will discuss several prescriptions in section 6. 

The HTL perturbation expansion generates ultraviolet divergences. In QCD pertur- 
bation theory, renormalizability constrains the ultraviolet divergences to have a form that 
can be cancelled by the counterterm Lagrangian A£qcd- We will demonstrate that the 
renormalization of HTLpt can be implemented by including a counterterm Lagrangian 
■A^HTL among the interaction terms in (2.3). There is no proof that the HTL perturbation 
expansion is renormalizable, so the general structure of the ultraviolet divergences is not 
known; however, it was shown in previous papers [46, 47] that it was possible to renor- 
malize the next-to-leading order HTLpt prediction for the pressure of QCD using only a 
vacuum counterterm, a Debye mass counterterm, and a fermion mass counterterm. By 
also including a coupling constant counterterm, this is possible at NNLO as well, as we 
will show in this paper. The necessary counterterms for the renormalization at NNLO as 
just discussed read 

A^o = + ^(^"«)) (1 - ^^'"^D , (2.5) 

Ami = {- ^^'u~J'^ "'^ + ^(^'"'0 " '^^"'^ ' ^^-^^ 
= (^-A-^as6 + 0(5^^?) (1 " , (2.7) 

SAas = _ 11^A-4^F ^2^2 ^ 0(^3^3) (2.8) 
127re 

where the vacuum and thermal mass counterterms were derived in ref. [46, 47], and the 
coupling constant counterterm is the standard one- loop running for QCD [56, 57]. 
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Dressed propagators 



Dressed vertex 




One-loop pure gauge contribution to gluon self-energy 



^^^= Gluon self-energy insertion 
^5^= Quark self-energy insertion 
^ Vertex insertion 




Figure 2. QCD diagrams contributing through NLO in HTLpt. The spiral lines are gluon propa- 
gators, the dotted lines are ghost propagators, and the solid lines are quark propagators. A circle 
with a n indicates a one-loop gluon self-energy insertion, a S indicates a one-loop quark self-energy 
insertion, and a F indicates a one-loop vertex insertion. A square with a 5 is shorthand for the pure 
gauge diagrams contributing to the one-loop gluon self-energy. All gluon and quark propagators 
and vertices shown are HTL-resummed propagators and vertices. The logic behind the diagram 
notation is as follows: diagrams consisting only of gauge propagators have g superscripts. Diagrams 
containing fermion propagators have / superscripts. The subscript indices are identical to those 
used in [49] (pure gauge QCD) and [45] (QED). We do not display the symmetry factors in the 
diagrams. 
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3 Diagrams for the thermodynamic potential 

Although HTLpt is defined in Minkowski space, it is much more convenient to carry out 
the necessary calculations in Euclidean space since in the paper we focus on static quanti- 
ties. In the imaginary-time formalism, Minkoswski energies have discrete imaginary values 
Pq = i2mTT for bosons and pq = i(2n + l)27rT for fermions, and integrals over Minkowski 
space are replaced by sum-integrals over Euclidean vectors (2n7rT, p) or (2(n -|- l)7rT, p), 
respectively, with n an integer. We will use the notation P = (Po,p) for Euclidean mo- 
menta. The magnitude of the spatial momentum will be denoted p = |p|. The inner 
product of two Euclidean vectors is P ■ Q = PqQo + P • q- The vector that specifies the 
thermal rest frame remains n = (1,0). Here we use dimensional regularization and the 
remaining three-dimensional integral is generalized to (i = 3 — 2e spatial dimensions. We 
define the dimensionally regularized sum-integral by 

•^■^ ^ ^ Po=2n7rT ^ ' 

•'^^^ ^ ^ Po = (2n+l)7rT ^ ^ 

where // is an arbitrary momentum scale. The factor (e"'^ /4irY is introduced so that, after 
minimal subtraction of the poles in e due to ultraviolet divergences, fi coincides with the 
renormalization scale of the MS renormalization scheme. 

The Feynman diagrams discussed in this section are gathered in figures 2 (notation 
key, LO, NLO) and 3 (NNLO). Using the same notation for the group theory factors as 
Arnold and Zhai [6], for QCD with Nc colors and Nf flavors of quarks, we have 

dA = N^-l, dF = N,Nf , CA = Nc, ^F = ^ , S2F= ^ Nf . 

(3.3) 

The earher work on NNLO pure-glue QCD [49] and NNLO QED [45] essentially treat 
special cases where the group symmetry factors are 

cIa = N^ — 1, dp = , CA = Nc , sp = , S2F = , (pure-glue QCD) 

(3.4) 

or 

dA = l, dF = Nf , CA = , SF = Nf , s2F = Nf , (QED) (3.5) 

and as such the analytical results from those papers will be highly useful in the present 
calculation as well. 

The thermodynamic potential at LO in HTLpt for QCD reads 

^^LO = dA [TL + Tf^] + dpTl^ + AoSo . (3.6) 

Here, J-'f^ is the contribution from the gluons and J-f■^^ the contribution from the ghost, and 
Aq^o is the leading-order vacuum counterterm. The expression of eq. (3.6) is presented in 
ref. [47], and the details can be found in section 3 therein. 
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The thermodynamic potential at NLO reads 



NLO 



2 d 2 "5 



77/ I 77/ 



2d 



(3.7) 



where AiiSq, Aimjj, and Aim^ are the terms of order 6 in the vacuum energy density and 
mass counterterms. Again the detailed expression of eq. (3.7) is presented in section 3 of 
ref. [47]. 

The NNLO HTLpt thermodynamic potential for QCD can be written as 



^^NNLO = ^^NLO + dACA 
+dACA 



773 I 779 I 779 _i_ 779 I 779 I 779 _i_ 779 



779 I 779 I 779 I 779 j_ 779 
•'ah ^ •'Si "I" •^ 3j "I" •^ 3k •^Sl 



+ dAS2F 



■pJ _|_ 77/ 
•^3a ^ -^3b 



+dASp 
+dASF 



He + H 



+ dACASF 



_jrJ i jrJ _|_ 77/ I 77. 



3m 



3n 



3o 



-CJ \ fJ \ -pJ \ fJ \ fJ I -fl 
-^3d ^ •^Se ^ •^Sf ^ •^Sg ^ •^Sk ^ •^Sl 



Q Q 



52 



(5m 



2 \2 
Dl 



LO 



2 \2 



[Aiml] 



+- 



1 



92 



2 [(am2)2 



(Aim2)%dA 



'^-4->^2a+2b+2c + '^-^•^2a+2b 



a. 



(3.8) 



where A2£^0; A.2m'jj, ^2171^ and A2as are the order 6"^ counterterms in the vacuum energy 
density, masses and coupling constant that can be read out from eqs. (2.5), (2.6), (2.7), 
and (2.8). 

The expressions for the bosonic diagrams J^ia~-^im presented in section 3 of ref. [49] , 
and the diagrams with fermions ^loT-^li ^'^^ expressed in section 3 of ref. [45]. The only 
NNLO diagrams specific to QCD, i.e. the non-Abelian diagrams involving quarks, are given 

by 



Hr. = \9'i Tr 

-3 T{PQR} 



r"(i? - P, R, P)SiP)T^{P - Q, P, Q)SiQ)T^Q - R, Q, R)SiR) 



3n 



3o 



,p/.^5^p _jlQ_pjl_ Q)A"^(P - R)A^^{Q - P)A^\R - Q) , 
n^''(P)A'^°(P)n°'^(P)A^^(P) , 

-Ig'i Tr \r"P{P,-P,Q,Q)S{Q)] A°'^(P)A^-(P)nf (P) , 
2 Tp{Q} L J 



(3.9) 
(3.10) 

(3.11) 
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where 



nf(p) 



1 2^ r^-'"'^(p,-p,Q, 



Q)A°''(Q) 



Tq 



Q)A'^^(Q)rT^(P, Q, -P - Q)A^\-P 



Q) 



Tq 




(3.12) 



nf (P) = Tr [r^(p, Q, g - p)5(Q)r'^(p, q, q - - p)] . (3.13) 



Thus Ii^^{P) is the one-loop gluon self-energy with HTL-resummed propagators and ver- 



4 Expansion in the mass parameters 

In refs. [46, 47], the NLO HTLpt pressure was reduced to scalar sum-integrals. It was clear 
that evaluating these scalar sum-integrals exactly was intractable and the sum-integrals 
were calculated approximately by expanding them in powers of nriD/T and mq/T following 
the method developed in ref . [30] . We will adopt the same strategy in this paper and carry 
out the expansion to high enough order to include all terms through order if m/j and 
ruq are taken to be of order g. 

The pressure can be divided into contributions from hard and soft momenta, which are 
the momenta proportional to the scales T and gT respectively. In the one-loop diagrams, 
the contributions are either hard (/i) or soft (s), while at the two-loop level, there are 
hard-hard (/i/i), hard-soft (/is), and soft-soft [ss) contributions. At three loops there are 
hard-hard-hard (hhh), hard-hard-soft (hhs), hard-soft-soft (hss), and soft-soft-soft (sss) 
contributions. As mentioned above, the (h), (s), (hh), (hs) and (ss) contributions have 
been calculated in refs. [46, 47], while most of the three-loop contributions can be read out 
from the earlier NNLO HTLpt work on pure-glue QCD [49] and QED [45], thus to keep 
the discussion compact, we only list the contribution that is specific to QCD here, i.e the 
term proportional to casf in eq. (3.8). 

If all momenta in the three loops are hard, we can obtain the m/j/T and niq/T ex- 
pansion by expanding in powers of and m^, and in this case to obtain the expansion 
through order g^, we can simply use bare propagators and vertices. The contributions from 
the three- loop diagrams were first calculated by Arnold and Zhai in ref. [5], and later by 



tices: 



ii^''{p) = CA^''{P) + SFnf{P). 



(3.14) 
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Braaten and Nieto in ref. [8] from which the [hhh) contribution to the casf term reads 



2 3a 



9\d- 



■f(hhh) 



f{hhh) 



3m+3n 



2{d-b)i 



1 



1 



^^"^ ^~^t{PQ}RP'Q'{P-QY{Q-RY 



864 



+2 

1 

e 



nPQR} P'QHP - Qm - RY{R - p? 

369 
250 



^1 2 + - 104 C^(-1) 2 C(-3) 
125 °^ ^25^-^ 25 C(-l) 25C(-3) 



47rr 



(4.1) 



Note that there is no {hhh) contribution from the HTL-specific diagram J^^^ [45]. 

In the {hhs) region, the momentum P is soft,^ while the momenta Q and R are 
always hard. The function that multiplies the static propagator Aj'(0, p), A/,(0,p) or 
Ax(0,p) can be expanded in powers of the momentum p. In the case of Ar(0, p), the 
resulting integrals over p have no scale and thus vanish in dimensional regularization. The 
integration measure scales like ml,, the static propagators A2,(0,p) and Ax(0,p) scale 
like 1/m^, and every power of p in the numerator scales like m£,.^ The {hhs) contribution 
to the CASp term reads 



■f(hhs) 



3a 



f(hhs) 



3m+3n+3o 



9^T{d 



p {p^ + m|,)2 
lid + 46) 




1 



Q{R} Q'^R'^ 

p2 



P 



p{p' + ml)^TQ{R} Q'R' 



p (p2 + 171%) 



3 



- H h 27£; 

e 21 ' 



Ar'^ 4 
38 



Tr 



R^\{R-Y)^/- 



KAttTJ \2mD J 



2e 



Sirrrip) 



(4.2) 



^The soft contribution arises fi'om the Po = term in the sum-integral. See appendix A of ref. [58] for 
details. 

^We refer the reader to e.g. ref. [46] for HTLpt Feynman rules, notations and conventions. 
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with the function Tp defined by 



where the angular brackets represent an average over c 



^^ = (^^2^)' (4-3) 



mni-e 



U{c)\ = ^ ':L. ' . / dc{l- c'rm ■ (4.4) 



For all of the diagrams that are infrared safe, the (hss) contribution is of order g'^m'jj, 
i.e. g^, and can be ignored. The infrared divergent diagrams, i.e. J'l^, contribute to the 
CASp term as follow 



-r-fihss) 



T{R} ^ Jpq 



2{p^ + 4m|,) 



5 The thermodynamic potential 

In this section we present the final renormalized QCD thermodynamic potential. The LO 
and NLO results were derived in refs. [46, 47], thus to keep the presentation compact, we 
only list the NNLO one, which corresponds to order in HTLpt. 
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Using the results listed in section 4, the renormalization contributions at order 5^ read 
d d d d 



dm 



D 



TO 



D 



+- 



+dA 



92 



LO 



2\2 



(Aim^) 



'^^-^2a+2b+2c + '^^-^2a+2b 



ideal ^ 



a, 
5 /I 



45/1 



, +2 bg- + 2 + 2-^7^ — f 
2 V e ^2 C(-l 



2 

45/1 



g +21og^ - 21ogmD + 2 ) "iD 



- + 2 + 21og I - 21og2 + 2||3^ ] 



+ 



315 / I 



235 /1 /i X'(-l) 149 132^ ^ 

- + 4 log — + 4^^^ — - H loe 2 

32 Ve 2 ^ C(-l) 47 47 ^ 



+ 



Spas 



+ - 



IT 

15 /I 



n 61 22, C'(-l) 

-+41og|-21ogmz) + — -ylog2 + 2^^^— ) m^j 

"25/1 u 12 , X'(-l) 



24 



_ + 41og c _ 21ogmi5 + ^ - 21og2 + '^^z^ 



(5.1) 



where Jideai is the pressure of a gas of massless spin-one bosons and tud, rhq and fi are 
dimensionless variables: 



-^ideal = 



TT 
'45' 



rhq 



2ttT ' 



27rr ' 
27rr ■ 



(5.2) 

(5.3) 
(5.4) 
(5.5) 



In order to keep the expression (5.1) compact, the pure-glue contributions have been rep- 
resented by Ar^y^ which is given by eq. (4.10) of ref. [49]. 

Adding the NNLO counterterms (5.1) to the contributions from the various NNLO 
diagrams we obtain the renormalized NNLO thermodynamic potential. We note that at 
NNLO all numerically determined subleading coefficients in e drop out and due to this the 
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final result is completely analytic. The resulting NNLO thermodynamic potential is 

1 



n 



NNLO 



ideal 



7 dp Spas 

4dA vr 



25 15 / [i 

Y + y"io + 15 (log I - ^ + + 21og2 ) ml. 



24 



+ 



Spas 



15 1 235 / u 144, 

log 

2 TTiD 16 V 2 47 



319 111 , 74C'(-1) 1 C'(-3) 

47^^+ 940 + 235 ^°S'-47CFr ^ 



315/ A 
+ ^(^°^2 
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log2 + 7£; + — mD + 90 



14 



47 C(-3) 

,2 



iriD 



+ 



Spas 
vr 



logrfiD 



5 1 
4 m£) 



+ 



25 
12 



, A 1 3 
^°^2 + 20 + 5^^ 



66, , , 4 C^(-1) 2CV3) 
^^°^'+5CM)"5CM) 



15 log 



1 



+ 7£; + 21og2 mD + 30 



niD 



+S2F 



as 



1 5 45 ' 

— (35 - 32 log 2) mfl 

64^ ^ ^ 2 



NNLO ' 



(5.6) 



where ^^^nlo NNLO pure-glue thermodynamic potential given by eq. (4.11) in 

ref. [49] We note that the coupling constant counterterm listed in eq. (2.8) coincides with 
the known one-loop running of QCD 



dfi 



{UcA - 4sp)g' 
247r2 



(5.7) 



In the next section we will present results as a function of g evaluated at the renormalization 
scale 27rT. 



6 Thermodynamic functions 

The mass parameters mpi and rUg in HTLpt are completely arbitrary. To complete a 
calculation, it is necessary to specify niD and as functions of g and T. As HTLpt is 
inspired by variational perturbation theory, the natural prescription for specifying the mass 
parameters is to use the variational solutions resulting from the gap equations 

d 

il.{T,as,mD,mq,ijL,6 = 1) = , (6.1) 

oniD 

d 

-^^VL{T, as,mD,mq, n,S = 1) = . (6.2) 

The LO gap equations are not well defined due to the lack of the coupling constant g in the 
LO thermodynamic potential, and in previous works on LO HTLpt one therefore chose the 
LO weak-coupling perturbative values of the Debye and thermal quark masses as physical 
estimates of the HTLpt mass parameters [42-44] . At NLO the variational prescription gives 
nontrivial solutions for the mass parameters, and the NLO results presented in [46, 47] use 
this prescription for both mp) and niq. For our NNLO results the variational prescription 
yields a complex Debye mass and niq = 0. The problem of a complex mass parameter was 
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Figure 4. Comparison of the real and imaginary valued contributions to the NNLO variational 
Debye mass, scaled by the LO perturbative Debye mass, with = 3 and Nf = 3. Note that in 
the small coupling limit the variational and LO perturbative Debye masses converge. 

encountered at NNLO in SPT massless 4>^ theory [29], in HTLpt QED [45], and in HTLpt 
pure-glue QCD [48, 49], so this seems to be a general issue with SPT/HTLpt at NNLO, 
which is not specific to QCD. Since the weak-coupling perturbative Debye mass receives 
contributions from the nonperturbative magnetic scale beyond LO, for the NNLO HTLpt 
pure-glue QCD results the Braaten and Nieto's (BN) mass parameter of three-dimensional 
electrostatic QCD (EQCD) was used as a substitute for the perturbative Debye mass, 
effectively discarding the nonperturbative contributions. The perturbative quark mass does 
not suffer from this issue, and as such there is bound to be some "mixing" of prescriptions 
for the mass parameters of the NNLO full QCD results unless one accepts the variational 
complex m£) value and sticks to the variational prescription for both parameters. In this 
section we will discuss and compare numerically the aforementioned prescriptions, before 
making a choice of prescriptions to use for comparison of our results to lattice data. 

6.1 Variational masses 

Solving the NNLO gap equations (6.1)-(6.2) yields a complex Debye mass and niq = 0. 
Figure 4 shows the real and imaginary parts of the complex NNLO Debye mass, scaled 
by the LO perturbative niD- While the imaginary contribution is non- negligible at inter- 
mediate coupling, and even grows to surpass the real contribution as approaches unity, 
it disappears in the small coupling limit, where the NNLO variational and LO perturba- 
tive results for niD overlap. Presumably the complex mass parameter is an artifact of the 
truncation after the m£,/T and rUq/T expansions, but this is difficult to confirm without 
extending the truncation to higher order. 

One strategy is then to throw away the imaginary part of the thermodynamic poten- 
tial to obtain thermodynamic functions that are real valued. Figure 5 shows the NNLO 
pressure found using this approach. For reference, figure 6 compares the real valued and 
the discarded imaginary contributions to the NNLO variational pressure over the same 
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Figure 5. The LO, NLO, and NNLO scaled variational pressures, with the imaginary contributions 
to the NNLO pressure discarded, Nc = 3 and Nf = 3. The shaded bands show the result of varying 
the renormalization scale /U by a factor of 2 around the central value n = 2'kT. 
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Figure 6. Comparison of the real and imaginary contributions to the NNLO variational pressure, 
with iVc = 3 and Nf = 3. The shaded bands show the result of varying the renormalization scale 
/J, by a factor of 2 around the central value n — 2ttT. 



temperature range. 



6.2 Perturbative masses 



At LO in the coupling constant g, the Debye mass is given by the static longitudinal gluon 
self-energy at zero three-momentum, mjj = 11^(0,0), i.e. 
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At NLO the weak-coupling perturbative Debye mass becomes logarithmically infrared di- 
vergent, reflecting the contribution from the nonperturbative magnetic scale g^T [59, 60]. 
The nonperturbative contribution was calculated in ref. [61] and reads 
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where ?n.mag is the nonperturbative magnetic mass. Using this mass prescription for the 
Debye mass would then require input from e.g. lattice simulations. For this reason we will 
not consider this prescription and instead will use the Braaten-Nieto prescription detailed 
in the next subsection. 

The NLO quark mass is infrared safe, and was calculated for A^c = 3 and Nf = 2 in 
ref. [62], 
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However, to our knowledge there has not been results available for general Aj, and in order 
to make the comparison with lattice data feasible, we use instead, whenever it applies, the 
Ac = 3 LO quark mass 
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6.3 BN mass 



The strategy of equating to the BN mass has earlier been applied to NNLO HTLpt 
Yang-Mills theory [48, 49]. Inspired by dimensional reduction, one equates the rriD mass 
parameter with the mass parameter of three-dimensional EQCD [7], i.e. = niE- This 
mass can be interpreted as the contribution to the Debye mass from the hard scale T and is 
well defined and gauge invariant order-by-order in perturbation theory. However, beyond 
NLO it will also depend on the factorization scale that separates the hard scale and the 
soft scale gT. In ref. [7] it was calculated to NLO, giving 
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Figure 7. Comparison of LO, NLO, and NNLO predictions for the scaled pressure for Nc = 3, 
Nf = 2 (left panel) and Nf = 3 (right panel), using BN and both LO perturbative niq and 
ruq = 0. 
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Figure 8. NNLO predictions for the scaled pressure with BN mo comparing LO perturbative rriq 
and niq = at renormalization scales \i = AttT and /i = ttT, for Nc = 3, Nf = 2 (left panel) and 
Nf — 3 (right panel). 



6.4 Choice of mass prescriptions 

To avoid dealing with imaginary contributions to the thermodynamic potential, we opt for 
the BN prescription for the Debye mass parameter, the same choice that was made for 
NNLO HTLpt pure-glue QCD [49]. 

As for niq, there is no problem with getting zero thermal quark mass from the varia- 
tional prescription: even if the m^-terms of Eq. (5.6) drop out, the quarks still contribute 
to the HTLpt thermodynamical potential through the terms proportional to dp, sp and 
S2F- On the other hand, when using the BN mass parameter one might want to use the 
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Figure 9. Comparison of LO, NLO, and NNLO predictions for the scaled pressure for Nf = 2 + 1 
(left panel) and Nf — 2 + 1 + 1 (right panel) with lattice data from Bazavov et al. [22] and Borsanyi 
et al. [23]. We use Nc — 3, three- loop running for as, n — 2ttT, and = 344 MeV. Shaded band 
shows the result of varying the renormalization scale /i by a factor of 2 around /i = 2t:T for the 
NNLO result. See main text for details. 

perturbative approach for the quark mass parameter, as using BN mass and perturbative 
thermal quark mass can be argued to be less of a "mixing" of prescriptions than BN mass 
and mq = 0. As it turns out, the final NNLO results are very insensitive to whether one 
chooses a perturbative mass prescription for niq, or uses the variational mass ruq = 0. In 
figure 7 the NNLO pressure with LO perturbative niq is virtually indistinguishable from 
the TTiq = result. The difference in dependence on the renormalization scale is also very 
small, though slightly in favor of the perturbative niq, see figure 8. However, convergence 
is improved with the choice niq = 0, where the pressure curves of figure 7 monotonically 
approaching the NNLO results rather than oscillating as one goes from LO to NLO to 
NNLO. We will therefore use rUq = in the following unless otherwise stated. 

6.5 Comparison to lattice data 

In figure 9 we show the normalized pressure for Nc = 3 and Nf = 2 + 1 (left panel), and 
Nc = 3 and Nf = 2 + 1 + 1 (right panel) as a function of T. The results at LO, NLO, 
and NNLO use the perturbative Debye mass given by eq. (6.7) as well as niq = 0. For 
the strong coupling constant as, we used three-loop running [11] with Aj^^ = 344 MeV 
which for Nf = 3 gives as(5GeV) = 0.2034 [12]. The central line is evaluated with the 
renormalization scale fi = 2ttT which is the value one expects from effective field theory 
calculations [7, 63] and the band represents a variation of /x by a factor of 2 around this 
value. 

The lattice data from the Wuppertal-Budapest collaboration use the stout action. 
Since their results show essentially no dependence on the lattice spacing (it is smaller than 
the statistical errors), they provide a continuum estimate by averaging the trace anomaly 
measured using their two smallest lattice spacings corresponding to N^- = 8 and Nt = 10 
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Figure 10. Comparison of NNLO predictions for the scaled trace anomaly with Nf = 2 + 1 (left 
panel) and Nf = 2 + 1 + 1 fermions (right panel) lattice data from Bazavov et al. [22] and Borsanyi 
et al. [23]. We use Nc — 3, three- loop running for as, n — 2ttT, and = 344 MeV. Shaded band 
shows the result of varying the renormalization scale /i by a factor of 2 around /i = 2TrT. See main 
text for details. 

[23], which were essentially on top of the Nj- = 6 measurement [64].^ Using standard 
lattice techniques, the continuum-estimated pressure is computed from an integral of the 
trace anomaly. The lattice data from the hotQCD collaboration are their Nr = 8 results 
using both the asqtad and p4 actions [22]. The hotQCD results have not been continuum 
extrapolated and the error bars correspond to only statistical errors and do not factor in 
the systematic error associated with the calculation which, for the pressure, is estimated by 
the hotQCD collaboration to be between 5 - 10%. We note that there are hotQCD results 
for physical light quark masses [65]; however, these are available only for temperatures 
below 260 MeV and the results are very close to the results shown in the figures so we do 
not include them here. 

As can be seen from figure 9 the successive HTLpt approximations represent an im- 
provement over the successive approximations coming from a weak-coupling expansion; 
however, as in the pure-glue case [48, 49], the NNLO result represents a significant correc- 
tion to the LO and NLO results. That being said, the NNLO HTLpt result agrees quite 
well with the available lattice data down to temperatures on the order or 2 Tc ~ 340 MeV 
for both Nf = 3 (left panel) and Nf = 4 (right panel). Below these temperatures the 
successive approximations give large corrections with the correction from NLO to NNLO 
reaching 100% near Tc. 

In figure 10, we show the NNLO approximation to the trace anomaly (interaction 
measure) normalized to as a function of T for Nc = 3 and Nf = 2 + 1 (left panel) 

*We note that the Wuppertal-Budapest group has published a few data points for the trace anomaly 
with Nt = 12 and within statistical error bars these are consistent with the published continuum estimated 
results. 
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and for Nc = 3 and Nf = 2 + 1 + 1 (right panel)."'' In the left panel we show data 
from both the Wuppertal-Budapest collaboration and the hotQCD collaboration taken 
from the same data sets displayed in figure 9 and described previously. In the case of the 
hotQCD results we note that the results for the trace anomaly using the p4 action show 
large lattice size affects at all temperatures shown and the asqtad results for the trace 
anomaly show large lattice size effects for T > 200 MeV. In the right panel we display 
a parameterization (solid blue curve) of the trace anomaly for Nf = 4 published by the 
Wuppertal-Budapest collaboration [23] since the individual data points were not published. 
In both the left and right panels we see very good agreement with the available lattice data 
down to temperatures on the order of T ~ 2Tc- 

6.6 Large Nf 

In the limit Nf ^ 1 while keeping g^Nf ~ 1, only ring diagrams contribute to the pressure 
in perturbation theory. Since is proportional to Nf in both QCD and QED, this indicates 
the equivalence of QED and QCD in the large- A^/ limit, with the large Nf effective coupling 
defined as ^eff = 9^/Nf/2 for QCD, and ^cff = e^/Nj for QED. The NNLO HTLpt results 
in this limit were discussed in the context of QED in ref. [68]. 

It is possible to solve for the 0{Nj) contribution to the pressure exactly [66], en- 
abling comparison between predictions from approximations such as perturbation theory 
or HTLpt, and exact numerical results. In figure 11, we plot the NLO and NNLO HTLpt 
predictions for the large- A^j pressure along with the numerical result of [66] as well as the 
perturbative g^g, g^g, and g^g predictions at ^ = e'^^nT obtained in [67]. For the pressure, 
going to NNLO in HTLpt extends the range of agreement with exact results compared to 
NLO, from g^s ^ 2 to ^2.8. The HTLpt large coupling behavior qualitatively matches 
that of the numerical results. 

7 Summary and outlook 

We have calculated the NNLO contributions to the thermodynamic functions of SU(A^c) 
Yang-Mills theory with Nf fermions using HTLpt. Using the BN mass for mo and niq = 
as the mass prescriptions, at NNLO we find that HTLpt predictions for the pressure and 
the trace anomaly are in agreement with lattice data for Nc = 3 and Nf G {3,4} down 
to T ~ 2 Tc. The failure of HTLpt to match lattice data at lower temperatures is to be 
expected since one is expanding around the trivial vacuum = and therefore neglects 
the approximate center symmetry Z(Nc), which becomes essential as one approaches the 
deconfinement transition [69-75]. In addition, it is also in line with expectations since 
below r ~ 2 — 3Tc a simple "electric" quasiparticle approximation breaks down due to 
nonperturbative chromomagnetic effects [59, 60]. 

^The Wuppertal-Budapest data were obtained using a physical charm mass at low temperatures; how- 
ever, we use massless quarks. The difference between massive and massless quarks is expected to be 
significant only for T < 414 MeV corresponding to the temperature where the lowest fermionic Matsubara 
mode equals the charm quark mass. 
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Figure 11. The predictions for the large- A^/ pressure of QCD between the numerical exact result 
from [66], NLO and NNLO HTLpt, and perturbative g^^ through gl^ [67] results at /it = e-^^nT. 



We find that when including quarks the agreement with lattice data is greatly improved 
as compared to the NNLO results of pure-glue QCD [48, 49]. Fermions are perturbative 
in the sense that they decouple in the dimensional-reduction step of effective field theory, 
so we expect that including contributions from quarks gives at least as good agreement 
with the lattice calculations as the pure-glue case. However, the exact reason for the better 
agreement between the HTLpt predictions and lattice calculations when including quarks 
is not clear to us. 

Just as for NNLO pure-glue QCD we found that the variational solution for the mass 
parameter m^i is complex and we therefore chose instead to use the Debye mass from 
EQCD for our vriD parameter. In the variational prescription the quark gap equation gives 
niq = 0, but even if the quarks do not get a thermal mass, the thermodynamical potential 
still gets HTLpt contributions from the inclusion of quarks in the calculation. Whether 
the complexity of the variational Debye mass parameter is due to the additional expansion 
in rriD/T and rUq/T is impossible to decide at this stage. The correction to the pressure 
going from NLO to NNLO is also rather large. It is unfortunate that the nonperturbative 
magnetic scale prevents going to N'^LO without supplementing the calculation with input 
from three-dimensional lattice calculations, as it would be interesting to see whether the 
complexity of the Debye mass parameter and the slow convergence persists. 

In closing, we emphasize that HTLpt provides a gauge invariant reorganization of 
perturbation theory for calculating static and dynamic quantities in thermal field theory. 
Given the good agreement with lattice data for thermodynamic quantities down to tem- 
peratures that are relevant for LHC, it would therefore be interesting and challenging to 
apply HTLpt to the calculation of dynamic quantities, especially transport coefficients, at 
these temperatures. 
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